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Abstract 

Discrete symmetries are commonplace in field theoretical models but pose a 
severe problem for cosmology since they lead to the formation of domain walls 
during spontaneous symmetry breaking in the early universe. However if one 
of the vacuua is favoured over the others, either energetically, or because of 
initial conditions, it will eventually come to dominate the universe. Using 
numerical methods, we study the evolution of the domain wall network for a 
variety of field configurations in two and three dimensions and quantify the 
rate at which the walls disappear. Good agreement is found with a recent 
analytic estimate of the termination of the scaling regime of the wall network. 
We conclude that there is no domain wall problem in the post-infiationary 
universe for a weakly coupled field which is not in thermal equilibrium. 
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I. INTRODUCTION 



It was first noted by Zel'dovich, Kobzarev and Okun [Q] that the restoration of sponta- 
neously broken discrete symmetries at high temperatures in the early universe poses severe 
problems for its subsequent evolution.^] If the manifold Ai of degenerate vacua of the theory 
is disconnected such that the homotopy group 7io{M) is non-trivial, then sheet-like topolog- 
ical defects — domain walls — form at the boundaries of the different degenerate vacuua 
during the symmetry breaking phase transition, due to the existence of a causal particle 
horizon in a decelerating universe The subsequent evolution of the wall network can be 
studied using the techniques of percolation theory and is such that the energy density 
in the walls eventually comes to dominate the total energy density. The consequence is an 
inflationary phase which suffers from a 'graceful exit' problem in that it creates an universe 
devoid of any matter [^|[. This can be avoided if the energy scale associated with the 
discrete symmetry breaking is low enough; however in order to avoid generating excessive 
anisotropy in the cosmic microwave background, the energy scale is further restricted to be 
smaller than ~ 1 MeV §]. Thus walls formed at any higher energy scale must be unstable 
and have decayed away long before the present epoch. 

Attempts to introduce instability in domain walls which are expected in field theoretical 
models have usually focussed on the role of non-renormalizable operators in the Lagrangian 
which may arise due to violation of global symmetries by Planck-scale gravitational effects 
0. These can 'tilt' the potential so as to favour one particular minimum which comes to 
dominate the universe because the pressure of the favoured vacuum eventually wins out over 
the restraining tension of the domain walls ||I|JlO|. An example from particle physics is the 
addition of a dimension-5 operator to the Lagrangian of the next-to-minimal supersymmet- 
ric standard model, in which the usual supersymmetric Higgs sector is supplemented by 
a gauge singlet superfield. Here Z3 domain walls are expected to be created in the Higgs 



fields during the electroweak symmetry breaking phase transition at T~mvi^ ]TT|. However 
the dissipation of the wall network in this case releases the contained energy as high energy 
interacting particles which can interefere adversely with cosmological processes such as pri- 



mordial nucleosynthesis |T2|, hence the walls must disappear well before this epoch [|T3|. In 
order to quantify this constraint, it is neccessary to study the rate at which the wall network 
dissipates for a given pressure difference between the different vacuua. 

Another way in which the domain wall network may be destabilized is through a suitable 
choice of initial conditions, viz. if the probability distribution has a bias for one vacuum 
over the others. This can be due to a prior non-equilibrium phase transition, for example 
cosmological inflation, which can displace weakly coupled light fields from the minima of 
their potential leading to a biased initial state [|T^]. It has been suggested that the 
resulting unstable domain walls in a light scalar field may be relevant for the formation of 
large-scale structure in the universe Recently, both an analytic |]T^ and a numerical 

181 study of such biased domain walls have been carried out. 



^It has been argued recently Q, following earlier work that symmetry restoration may not 
neccessarily occur, particularly if the scalar field has no gauge interactions. Whether this is still 
true when non-perturbative effects are taken into account is currently under investigation [Q]. 
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In this paper we make a systematic investigation of domain wall network evolution for 
both cases, viz. 'pressure' (§111 (A)) and 'bias' (§111 (B)) and and compare our results 
with simple analytic estimates of the rate at which the walls disappear, as well as with 
previous numerical work |[l7| , p!8[| . We also study (§111 (C)) the effects of biased initial field 
configurations resulting from a period of primordial inflation. We conclude (§IV) with a 
discussion of the physical implications of our results. 



II. NUMERICAL TECHNIQUES 

The generic potential which exhibits the problem under discussion is 

V(,)^V.[^_-l)\ (1) 

This has two degenerate vacuua, (j) = ±0o separated by a potential barrier Vq, so that the 
spontaneous breaking of the Z2 symmetry when the universe cools below the temperature 
T ~ Vo will result in the formation of domain walls. We wish to simulate the evolution of 
the resulting wall network by solving the field equations on a specified lattice. The basic 
problem in such numerical studies is that there are two very different length scales involved, 
viz. the wall width. 

Wo ~ , (2) 

and the size of the simulation box. The first is constant in physical coordinates whilst the 
second must be large relative to a typical domain size, which scales with the expansion of the 
Universe. However the wall thickness is in general much smaller than the size of of the wall 
network so one can sensibly assume that the walls behave like two-dimensional relativistic 
membranes, given that their internal structure is expected to be of negligible importance 
for the evolution. Two routes have so far been pursued in applying this idea to computer 
simulations of domain wall djTiamics. 



Kawano |]T9| considered an effective action obtained by expanding in wq/R, where R is 
the radius of curvature of the wall, and retaining only the zeroth order term. The resultant 
Nambu-type action yields the evolution equation 



R + 3^R{1 - R') = _^iL_^ , (3) 
a R 

where a is the cosmological scale-factor of the Robertson- Walker metric for an Einstein- 
DeSitter universe, g^i, = diag[— 1, a(t)^, a(t)^, a(t)^]. Unfortunately this otherwise attractive 
approach leads to severe numerical instabilities (of the type discussed in a slightly different 
context 1^), hence cannot be fruitfully pursued. 

The second approach, which we use in the present work is due to Press, Ryden and 
Spergel |22|,^. Consider the classical equation of motion for the Higgs field in an expanding 



background: 



d'^cj) ^ ^dlnal d(f) 

dr]"^ d\nr]r] dr] 90 ' 
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where t] is the conformal time (d?7 = dt/a{t)) which measures the comovm^^ distance traversed 
by hght since the big bang. Here the numerical problem is manifest in the term on the rhs 
which makes the potential barrier appear higher as time goes on, resulting in a wall solution 
which grows increasingly narrow (in comoving cordinates) with time. In attempting to 
evolve an equation of this type directly one would find that the wall solutions in the field 
quickly become so narrow as to be unresolvable. The solution is to generalize the equation 
of motion to: 

(9^0^ din a 190 is'^^ /g\ 

dr]"^ dlni] 7] dr] dcj) ^ 

and then set /? = in order to freeze the wall size in comoving coordinates. We also set 



a = 3 to ensure momentum conservation since this requires that we have [22 



« + f = 3. (6) 

Press et al. have discussed in some detail the justification for this approximation and 
demonstrated that it has a negligible effect on the evolution by performing simulations 
comparing a = (3 = 2 and a = 3, /5 = results over a limited range of rj. We perform a 
test of a complementary nature as described below which supports their contention that the 
numerical method is robust and unlikely to give spurious results. 

Following ref. |22| we implement this procedure using the standard finite differenceing 



scheme embodied in the expressions 

^ - \<-T^, > (7) 

(V^0)jj,fc = 0i+l,j,fc + 4>i-l,j,k + 4>i,j+l,k + 4>i,j-l,k + 4>i,j,k+l + 4>i,j,k-l — Q4'i,j,k , (8) 

;n+l/2 _ (l-^)C,fc°-'+A^(V^'A5fc-aW805J /Q^ 

<t^-^k = <t^i,k + A#i;r-'- (10) 



We also use the same algorithm |2^ for measuring the total wall area in the field. Our 
simulations were run in both two and three dimensions, on a 128x128x128 grid {D = 3) 
and a 1024x1024 grid {D = 2). 



III. WALL NETWORK EVOLUTION 

The general evolution of the domain wall network is well known ; the most interesting 
phase is that of 'Kibble scaling', in which the walls have a correlation length in comoving 
units given by 

^-vrj, (11) 

where v is the velocity of wall propagation. The simplest variable which tracks the evolution 
is then the comoving area of the wall network per unit volume, (A/V), which scales as 

(A/V) ^ (A/\/)o| ^ , ^° , , (12) 
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while the physical energy density contained in walls behaves as a ^ times this. We find that 
there is indeed a scaling region with 

/ 0.95 ±0.08 forD = 2, 
^^1^)^^ ' 0.92 ±0.08 forD = 3, ^^^^ 



in agreement both with Eq.(0) and other simulations |2^JT^. (The errors quoted here are 
purely statistical). Hence we can calculate f , as defined in Eq. ([TT|) , to be 

r0.7±0.1 forZ^ = 2, 

\ 0.5 ± 0.1 for = 3 , ^ ' 

(Note that the velocity which was found previously to be ~ 0.4 in three dimensions [32 



was 



somewhat differently defined.) We shall henceforth simplify our algebra by normalising the 
values of the comoving energy density and scale factor to be unity at the initial time of our 
simulation, here denoted with the subscript 0. The evolution of the network thus proceeds 
as follows. For the initial few timesteps, with ^ 10, the correlation length is less than the 



wall thickness, before the network enters the scaling regime of equation (12) which ends at 
?7 ~ 100 when it becomes comparable to the size of the lattice. 

To verify the suitability of the numerical method used we ran a suite of simulations for 
D = ?) where /3 is increased in steps of 0.1 from zero and a is adjusted according to Eq.(|D. In 
order for the domain wall to be resolvable, it should occupy at least one (preferably several) 
lattice spacing(s); this corresponds to a limiting value of /3 ^ 0.66 in our simulations. In 
figure |I| we show that the exponent v in Eq.(|T3p does not change significantly as /3 is increased 
up to this value.0 Overall the gross features of the evolution appear to be insensitive to the 
numerical approximation used; in particular the onset and termination of the scaling regime 
are not dependent on /3. This adds further weight to the conclusion of Press et al. that there 
is no difference between the dynamics of 'real' walls which have constant physical thickness, 
and the 'artificial' walls in these simulations which have constant comoving thickness. 



A. Pressure 



A well known way |TO|J21] to evade the causal limit on the disappearance of the domain 
walls is to include a pressure term in the potential, viz. 



V^(0) = Vo 



00 



(15) 



The dynamics is now expected to be dominated by two competing forces, the surface tension 
cr/R where i? ^ aC, is the radius of curvature of the wall structure, and the pressure which 



^There appears to be a small systematic trend of decreasing u with increasing /3. Although 
not statistically significant, this may just reflect the increasing efficiency of the area measuring 
algorithm at picking up small bubbles at late times, as P is increased. 
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is equal to the differences in energy density of the two minima of the potential, SV. Then 
we expect simple scaling behaviour (12) until some conformal time % at which the pressure 



becomes comparable to the surface tension, and the wall network disappears exponentially 
fast. Now rjc is given by 

a , , = SV . (16) 

where the subscript indicates the value at the begining of the evolution. In the present 
case we are scaling the fields with time so as to remove factors of a and set tjq = 1 so this 
equation will reduce simply to 

- = 2/i . (17) 

Vc 

We display this behaviour in figure H, where the comoving area density is plotted against 
conformal time for several choices of e. As expected, there is an exponential fall-off from the 
scaling regime at some value of rjc which decreases with increasing /i. The 'bounces' in the 
curves occur as the bubbles of the disfavoured vacuum collapse and radiate away the energy 



contained in the walls in the form of Goldstone bosons; as noted earlier [|25| a large fraction 
of the energy is lost in a few bounces. The relation between /i and rj^ is shown in figures |^ 
and ^ which show excellent agreement with the expected behaviour ( [T7| ) to within a factor 
of 2. (Here 7]^ is defined as the value of rj at which the product ri{A/V) has fallen off by 
some factorm taken here to be 10 and 100 for illustration.) The actual rate of the decay of 
the wall network after the exponential decay has set in is hard to measure. We find here 
that assuming a behaviour of form {A/V)c ~ r;"^ exp[— ^(/i?])"] with k some constant, n is 
^ 2 ± 1 (D = 2) and ^ 3 ± 1 (0 = 3). 

The implication of this is that we can regard the walls as simply disappearing as a 
consequence of the difference in energy between the two minima of the potential 5V at a 
correlation length Rc = a/6V (in physical, not comoving, units), where again the constant 
of proportionality is around unity. Hence, if we have walls forming at the electroweak scale 



and we require that they disappear before nucleosynthesis ||12|, we find 

— > -4 10-'" GeV , (18) 

a 1010 cm ' ^ ^ 

in good agreement with previous estimates made from physical arguments as to the time of 
pressure domination M] and 2-dimensional thin wall simulations ||2H . 



B. Bias 



Next we consider the possibility of a deviation from scaling behaviour due to a biased 
initial probability distribution. For concreteness we consider the Z2 case where there are 
only two distinct vacuua and we generate the initial configuration with a probability 



0.5 + e 



(19) 



that each initial domain is in the + phase, where e > 0. A similar exercise has been 



performed recently by Coulson, Lalak and Ovrut [T^ who find that the domain wall network 
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then evolves much more rapidly than in the usual case of Eq. (0) , and the favoured domain 
rapidly dominates the universe. Such a scenario is interesting because a bias in a light scalar 
field which is not in thermal equilibrium can be generated naturally by a previous epoch 
of inflation. Other possible ways in which a scalar field is more likely to find itself in one 
minimum than the other are if the symmetry is approximate, so that the two minima are not 
exactly degenerate, or if the symmetry breaking occurs through some intermediate phase 
which allows one minimum to be preferentially populated. 

The evolution is shown in figure ^ where we plot comoving area density against conformal 
time for several choices of e, ranging upto 0.03. Typically, the behaviour is similar to the 
usual case of no bias, with an exponential decay of the comoving energy density starting 
at some critical timescale rjc- (For rj ^ 10 the behaviour is unphysical since ^ is less than 
the domain wall thickness (here scaled to be 5), while for t] ^ 100, ^ is approaching the 
physical size of the box, and hence the behaviour is again unphysical.) Notice that the wall 
network dissipates rapidly as e is increased above zero, as was also found in ref. |T^. (Again 
we see the 'bounces' associated with the radiating away of the energy contained in the wall 
network.) 



Hindmarsh |jT^ has recently made an analytic study of biased domain walls. He considers 
the evolution of the domain walls in terms of a scalar field constructed so that its zeros 
correspond to the locations of the domain walls, and whose evolution can be calculated in 
terms of Gaussian average field configurations; the expectation for the wall surface area 
density in D dimensions is 

(A/V) ^^exp[-Ke^r]''), (20) 

where k is some constant which should, in our units, be of order unity. The above result 
was obtained using rather sophisticated techniques but we can derive it from a much simpler 
(but perhaps less trustworthy) counting argument. We consider the wall evolution to proceed 
in such a way that the universe naturally divides itself up into domains of size C,—vt] (in 
comoving coordinates). Such domains have been causally connected, and so have had time 
enough to organize themselves to be either all + or all — . Hence we can see that the 
comoving area density in the absence of bias will behave as [A/V) ~ ~ [vr])^^ as 
expected. We can now examine the behaviour of the comoving area density in the more 
complicated situation when there is a bias. Now we expect a domain of conformal size ^ to 
become a + domain if most of its subdomains are +, and a — domain otherwise. (We 
have normalized so that the size of the initial domains is unity.) If the probability of each 
subdomain entering the positive minimum is p+, the probability for the whole domain to be 
— is given by a sum over values of a Poisson distribution which can, in any interesting case, 
be taken to be a Gaussian integral, viz. 

^(domain of sites, mostly -) = erf ~^/^ \ ^ eTf{V26y/N) , (21) 

\ (Tn J 

where N = {0.5+e)N is the mean number of + sites and = Vn is the standard deviation. 
Here the Gaussian integral erf(x) = dte~* /^/-\/27r can be adequately approximated for 
large x by 
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erf(a;) = ^^e'^'/l (22) 

We see that for a box with ^ sites, the probabihty given in Eq.dll) is just erf (725^^/2) 
and hence we expect exactly the exponential behaviour shown earlier in Eg. (|20|) . The extra 
factor of {^^^'^)~^ is cancelled by a combinatoric factor from the many different ways in 
which we can lay down our domains upon the overall network. This uses the fact that a 
cluster containing sites will typically scale as N^^'^, using the results of percolation theory 
as recently applied to similar problems in cosmology [P^. 



To verify this behaviour, we define rj^ to be the conformal time at which the product 
rj{A/V) has decreased by some substantial factor, which we choose to be 10 and 100 for 
illustration. The values of i]'^^"^ are plotted against e for both 2 and 3 dimensions in 
figures ^ and ^ and we find good agreement with the theoretical prediction of a straight 
line. Note that for both cases, we have cut off the figures at ~ 10 since before this 
conformal time the domain wall network is not yet fully formed. We also expect unreliable 
results when t] reaches 100 or 1000 in three and two dimensions respectively, since then the 
correlation length is approaching the size of our whole simulation. While we can measure 
the relationship between rjc and the bias e, it is again harder to measure the actual rate of 
the exponential fall off, i.e. the exponent of rj in the exponential. We find this exponent to 
be ~ 1 — 2 {D = 2) and ^ 2 — 3 (-D = 3), in acceptable agreement with theory. 

Therefore we can, for cosmological purposes, regard the wall evolution as being in the 
Kibble scaling regime for early times, and modelled subsequently by an exponential fall off 
at some conformal time rjc, which will lead to an almost immediate collapse of the wall 
network. The required value of e is readily calculable in three dimensions from 

'3/2 ,m , -3/2 

, (23) 




where rjc (Tc) is the scale factor (temperature) at the time of disappearance and rjo (To) at 
the time of domain formation. (Note that here we put k ~ 1 from our simulation result.) 
For example, in the event of wall formation at the weak scale of order 100 GeV, with the 
requirement that the walls disappear before the nucleosynthesis era which starts at ~ 1 MeV 
|l^, we find that r]c/f]o ~ 10^, and hence that the bias in the initial probability distribution 



must be of order 10 ^ or greater. However if the walls form subsequent to reheating following 
inflation, at an energy scale of say 10^ GeV, the bias need only be about 10"^'*. 



C. Realistic Field Configurations 



It is usually assumed that following cosmological discrete symmetry breaking in a 
scalar field, the field values will be uncorrelated on scales larger than the causal (particle) 
horizon, so the degenerate vacuua will be populated with equally probability. However, a 
prior period of cosmological inflation can result in correlations on (apparently) super-horizon 
scales if the field is sufficiently weakly coupled so as not to be in thermal equilibrium. This is 
because quantum fluctuations during inflation (with Hubble parameter H) will induce long- 
wavelength fluctuations in all scalar fields with mass m <^ H on spatial scales ^ 



8 



jl^ . For a minimally coupled field with vanishing potential (e.g. a (pseudo-) goldstone bo- 
son), this leads to the formation of a classical inhomogeneous field, with gaussian probability 



distribution |E8 



V{(p) oc exp 



^ h — 



IP 



dln/c . 



(24) 



where (pk is the mean value of the field averaged over a domain of size and kc is an 
effective ultraviolet cutoff imposed by physical considerations (dependent on the nature 
of the field). Since infiation blows up the scale factor exponentially fast, this provides a 
natural mechanism for bias since the global ensemble average of may not be realized even 
over regions as big as our present universe. During the post-infiationary phase, the above 
distribution (p^) averaged over a chosen scale thus becomes skewed [|15| . We parametrize this 
for the Z2 case by drawing probabilities for populating the two vacuua from the distribution 



P(0) 



27Ta 



(0.5 -£)e 



-(.^+1)2/2^2 



+ (0.5 + e)e 



(25) 



where e is the effective bias over the chosen scale. (Presumably its value can be calculated 
in principle if the field and infiationary parameters are specified.) In figures |^ and ^ we 
show the comoving area density plotted against conformal time for several choices of a with 
e = 0, 0.005. Surprisingly, the exponential fall-off from the scaling regime appears to be 
quite insensitive to the width of the probability distribution. This supports the suggestion 
||T5| that this is an efficient way to eliminate the domain wall network. 



IV. CONCLUSIONS 

We have considered two possible ways in which a cosmological domain wall network 
can be made unstable. The first of these is that of pressure, i.e. a small breaking of the 
degeneracy between the minima. The physical implications of this are well understood and 
we have simply confirmed numerically the usual argument ||10| that when the typical domain 
size R has grown to a critical value -Rc, the wall energy density decays exponentially fast, 
so the network disappears essentially instantaneously. The value of Rc is 

Rc = ^ , (26) 
dp 

where a is the wall surface energy density and 6p is the difference in energy density between 
the two minima. (There may be an additional effect due to the bias induced at the phase 
transition by the energy difference between the two minima.) 

The second mechanism is that of bias, where the initial field configuration is not sym- 
metric between the various minima of the potential. Here we confirm previous analytical 



p!7| and numerical ||T8[ studies which predict an exponential decay in the energy density 
with a characteristic conformal time rjc. This is an attractive solution to the domain wall 
problem if we can find a natural mechanism to bias the initial distribution, particularly since 
we find the evolution to be insensitive to the width of the distribution. Since such a bias 
must appear on (apparently) causally disconnected scales, it can only occur if either the 
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minima are genuinely inequivalent (and so have non-degenerate energy densities) or else if 
we are in a post-inflationary phase in which the field has correlations on super-horizon scales. 
This will allow the universe to more efficiently organise itself into one preferred minimum 
everywhere, thus enabling the wall network to decay away. Now an initial field configuration 
with correlations on all wavelengths can be characterized by an effective bias s{R) averaged 
over a box of size R, as high frequency modes will average out to zero. We have seen that 
a region of size R will decay away exponentially after a conformal time i]c{R) ~ e{R)~^f'^. 
For a radiation dominated universe, the time t{R) at which the domain will decay is then 

t{R) ~ e{R)-^ , (27) 

which must grow more slowly than R to ensure that the ordinary causal scaling bound is ex- 
ceeded. Hence for bias to provide a realistic means of eliminating wall networks through an 
initially non-equilibrium field distribution, we must have that s{R) falls off with R more 
slowly than R~^/^. Since the spectrum of fluctuations from inflation is (nearly) scale- 
invariant, we may expect that the bias £{R) does stay approximately constant with i?, 
and hence that this mechanism provides a way of eliminating domain walls in flelds which 
are sufficiently weakly coupled so as not to be in thermal equilibrium. 
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FIG. 1. A test of the numerical approximation used. The scaling index u of the wall network 
is seen not to vary significantly as the parameter /3 in the equation of motion is varied. 
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FIG. 3. ri~ against pressure /x in 2 dimensions, where rjc is the conformal time at which the 
product of r) and the comoving area density has fallen by a factor 10 (circles) or 100 (squares). 
The line connects the mean value of r)c over all the runs. 
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FIG. 4. r]~^ against pressure in 3 dimensions, where rjc is the conformal time at which the 
product of r] and the comoving area density has fallen by a factor 10 (circles) or 100 (squares). 
The line connects the mean value of r)c over all the runs. 
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FIG. 5. Comoving area against conformal time in 3 dimensions with the bias e in the range 
- 0.03. 




FIG. 6. 77~^ against bias e in 2 dimensions, where r/c is the conformal time at which the product 
of r] and the comoving area density has fallen by a factor 10 (circles) or 100 (squares). The line 
connects the mean value of 77c over all the runs. 
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FIG. 7. r]~^ against bias e in 3 dimensions, where rjc is the conformal time at which the product 
of r) and the comoving area density has fallen by a factor 10 (circles) or 100 (squares). The line 
connects the mean value of r)c over all the runs. 




FIG. 8. Comoving area against conformal time for post-inflationary field configurations, 
parametrized by the width a = 0, 0.5, 1, 1.5, 2 with no bias. 
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FIG. 9. Comoving area against conformal time for post-inflationary fleld configurations, 
parametrized by the width (7 = 0, 0.5, 1, 1.5, 2 with bias e = 0.005. 
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